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Abstract Granular materials transmit stress via a net¬ 
work of force chains. Despite the importance of these 
chains to characterizing the stress state and dynamics 
of the system, there is no common framework for quan¬ 
tifying their their properties. Recently, attention has 
turned to the tools of network science as a promising 
route to such a description. In this paper, we apply com¬ 
munity detection techniques to numerically-generated 
packings of spheres over a range of interparticle friction 
coefficients and confining pressures. In order to extract 
chain-like features, we use a modularity maximization 
with a recently-developed geographical null model [T], 
and optimize the technique to detect branched struc¬ 
tures by minimizing the normalized convex hull of the 
detected communities. We characterize the force chain 
communities by their size (number of particles), net¬ 
work strength (internal forces), and normalized convex 
hull ratio (sparseness). We find the that the first two ex¬ 
hibit an approximately linear correlation and are there¬ 
fore largely redundant. For both pressure P and inter¬ 
particle friction /r, we observe crossovers in behavior. 
For /i < 0.1, the packings exhibit more sensitivity to 
pressure. In addition, we identify a crossover pressure 
where the frictional dependence switches from having 
more large/strong communities at low /i vs. high ja. We 
explain these phenomena by comparison to the spatial 
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distribution of communities along the vertical axis of 
the system. These results provide new tools for con¬ 
sidering the mesoscale structure of a granular system 
and pave the way for reduced descriptions based on the 
force chain structure. 

Keywords force chains • network community struc¬ 
ture • numerical simulations • friction 


1 Introduction 

For more than half a century, it has been common to 
visualize the heterogeneous force transmission in gran¬ 
ular materials misiiii: these patterns have come to be 
known as force chains. It has been less clear, however, 
how best to provide a mesoscale description of this net¬ 
work of interparticle contacts. A better understanding 
of the important length scales over which intermedi¬ 
ate structures are present would provide new routes to 
connect particle-scale properties to bulk properties. In 
this paper, we take the growing field of network sci¬ 
ence as our inspiration [5]. Network techniques can be 
applied to such varied systems as social networks, neu¬ 
ral systems, or airline route maps: anything which can 
be reduced to a network of nodes and the edges (links) 
that connect them. In the sphere packings studied here, 
the network is composed of a set of particles (nodes) 
and the normal contact force between adjacent parti¬ 
cles (edges). 

The use of network science techniques has attracted 
significant attention in the past few years, particularly 
as a way to extract the “backbone” of the most impor¬ 
tant particles in the system, and to follow the evolution 
of those networks under loading. By considering the en¬ 
tire network miniiTiiiii] it is possible extract statistics 
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about the degree of connectivity among the particles, 
and how that influences the bulk response. 

One approach has been to define the main network 
based on a set of rules about how particles in a force 
chain should be connected. For example, Peters et al. 
m and Zhang et al. m define the backbone by set¬ 
ting a threshold value for contact forces and the angles 
between particles, while Kondic et al. [12] use a topo¬ 
logical invariant called the zeroth Betti number [13] to 
characterize the size of connected clusters. A disadvan¬ 
tage of these approaches is that thresholding strictly 
removes the weak interparticle contacts from consider¬ 
ation, even though of these forces are thought to play 
an important role in providing lateral stability [T4] . 
To address this issue, it is possible to use community- 
detection techniques [THUS] that allow for optimized 
partitioning into clusters without a hard threshold. For 
example, Navakas et al. [m HH] and Bassett et al. [3] 
identified force clusters which have stronger interparti¬ 
cle forces within each cluster than between them. How¬ 
ever, the force chain communities detected in this way 
take the form of compact domains rather than branched 
networks. Therefore, they seem unlikely to be the cor¬ 
rect mesoscale units from which to explain sound prop¬ 
agation [19] (observed to be along force chains) or the 
changing network of contacts under flow [20] (they form 
a giant component of broken links). 

Recently, Bassett et al. [T] recognized that community- 
detection algorithms which depend on a random null 
model miss an important aspect of granular materials: 
that grains are geographically constrained to be con¬ 
nected only to their neighbors. By working with a null 
model that respects these geographic constraints, the 
detected communities take the form of the expected 
branched structures. This geographical null models when 
applied to either simulated (frictionless) or laboratory 
(frictional) packings of disks, was able to distinguish 
the different force chain morphologies of the two dis¬ 
tinct datasets. 

In this paper, we adapt this technique to perform a 
similar community detection in simulations of 3D gran¬ 
ular packings, and examine how these communities sys¬ 
tematically change as a function of confining pressure 
and interparticle friction. By using simulations, we can 
test the methods in a controlled environment where it is 
possible to generate many independent realizations. In¬ 
creasingly, such data is becoming available in 3D gran¬ 
ular experiments [21] [22] [23], in the form of normal 
contact forces measured from the macroscopic defor¬ 
mations of soft particles. 

The community detection method consists of two 
steps, modularity maximization to partition the net¬ 
work into clusters, and selecting a resolution parameter 


P=0 P=10'® P=3xl0'^ 



Fig. 1 (A) Sample particle configurations generated using 
LAMMPS, with small particles (diameter d) shown in blue 
and large particles (lAd) shown in yellow. In subsequent 
steps, we apply pressure from above using a flat slab. (B) 
Corresponding normal force networks, with bar thickness pro¬ 
portional to the normal force magnitude. 

which controls the total number and shape of these clus¬ 
ters. To perform the maximization, we use the same ge¬ 
ographic null model as in Bassett et al. [T], allowing us 
to incorporate contact information. In selecting a reso¬ 
lution parameter, we found that the previous technique 
[1] was inadequate for 3D systems. Therefore, we de¬ 
veloped a new figure of merit to quantify the degree of 
branching within the communities, based on the convex 
hull of the constituent particles. 

When the process is complete, each packing is par¬ 
titioned into a set of branched communities (with many 
interstitial communities consisting of only a few weak- 
force particles). We characterize the ensemble of com¬ 
munities by their size, strength, and degree of branching 
as a function of both interparticle friction and pressure. 
Both friction and pressure influence the network prop¬ 
erties of force chains, and we observe crossovers as a 
function of both parameters. 

2 Simulation Methods 

We perform our numerical simulations using the dis¬ 
crete element model LAMMPS (Large-scale Atomic/ 
Molecular Massively Parallel Simulator) [24] maintained 
by Sandia National Labs. This open-source software is 
based on a fast parallel algorithm [25] for molecular 
dynamics. Our simulations contain N = 3000 bidis- 
perse spheres poured from above, half of diameter d 
and half of diameter lAd. The simulation cell has lat¬ 
eral dimensions Ibd x Ibd (periodic boundary condi¬ 
tions in both directions) and height 26d (open at the 
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top, closed at the bottom) with gravity acting down¬ 
wards. To model a hard, frictional granular material, 
we use a Hertzian contact model with a normal elastic 
constant = 2 x 10 ^ and a tangential elastic constant 
Kt = |i^n- These values, when re-dimensionalized, ap¬ 
proximately correspond to ruby spheres [26] of centi- 
metric size. We vary the interparticle friction coefficient 
over seven values /i = 0, 0.01, 0.03, 0.1, 0.3,1, 3 to exam¬ 
ine the dependence of our results on interparticle fric¬ 
tion. 

We mimic experimental protocols in which particles 
are poured into a box from above, and then compressed 
via a uniform pressure. Our numerical pouring method, 
adapted from Silbert et al. m , mimics pouring par¬ 
ticles through a sieve (to prevent the formation of a 
conical heap) by generating particles at random posi¬ 
tions at the top of the simulation volume and allowing 
them to fall downward under the force of gravity. We 
select a very low packing fraction (j)i = 0.005 for this 
insertion region, so that the resulting mean coordina¬ 
tion number c = 5.5 is independent of the choice of <pi 
[28] . After inserting all particles into the container, we 
allow the kinetic energy to dissipate until it is less than 
l^-^mgd [27]. 

For each such initial packing, we apply pressure by 
generating a massless slab (size 15d x 15d x 2 (i, FCC 
lattice, as shown in Fig. at the top of the packing. 
From this initial state (compressed by gravity, but no 
additional pressure) we apply an increasing series pres¬ 
sures (P = 10-^3 X 10-^10-^3 X 10-^10-^, mea¬ 
sured in units of K^) to the upper surface of the result¬ 
ing packing. The magnitude of the non-dimensionalized 
confining pressure is just below those reported in recent 
experiments on softer particles [29]. After the kinetic 
energy of the system has again dissipated, we record po¬ 
sition and force measurements before iterating through 
all 6 pressure values. At each step, we record all particle 
positions and interparticle forces; sample normal force 
networks are shown in Fig. 

For each of the seven values of /i, we perform 20 in¬ 
dependent simulations starting from different random 
initial conditions. We use a bootstrap-like process (sam¬ 
pling with replacement) to confirm that this is sufficient 
for reliable statistics. In a few places, noted within the 
text, the fluctuations were large enough that this crite¬ 
rion was not satisfied. In all of our analyses, we consider 
only the normal component of the interparticle forces, 
as is currently measured in experiments [ 2 TI [22l [23]. 
This simplification also allows us to directly compare 
both frictionless and frictional packings. 


3 Community Detection 


Our goal is to partition the granular packing into com¬ 
munities of particles which have high interparticle forces 
internally (locally stiff) and low interparticle forces in 
their connections to other communities. We build on 
the work of Bassett et al. HE], which utilizes the open 
source network analysis tool GenLouvain (Version 2.0) 
from Net Wiki m to implement the modularity maxi¬ 
mization method m of community detection. 

The network analysis begins from a representation 
of the normal force network (Fig. Ef) as an V X V 
weighted adjacency matrix W. Each element Wij is 
zero for particles not in contact, and fij/f for all non¬ 
zero interparticle forces (scaled by the mean normal 
force / for the whole packing). The modularity Q of a 
network is a scalar value calculated from 


Q — [^ij 7^ij] ^j) (1) 

hj 


where 7 is a resolution parameter, Pij is the expected 
weight of an edge due to a specific null model, Ci and 
Cj are the (numbered) community assignments for par¬ 
ticle i and j, and 6 is the Kronecker delta function. If 
particles i,jf are assigned to in the same community, 
then 6{ci,Cj) = 1, otherwise 6{ci,Cj) = 0. The opti¬ 
mization process adjusts the community assignments 
for fixed 7 and fixed null model. As developed in Bassett 
et al. [ 1 ], we utilize a physically-motivated geographic 
null model in which particles connect to a community 
through their direct neighbors: 



Wij ^ 0 , 

Wij = 0. 


( 2 ) 


We chose not to use the more common Newman-Girvan 
null model [31] which allows for arbitrary connections 
between particles mm because in 2D packings of disks, 
the geographic null model has been shown to success¬ 
fully generate communities with chain-like morpholo¬ 
gies [ 1 ] (rather than compact domains). 

Because modularity maximization (finding the largest 
value of Q) is an NP-hard problem, the published meth¬ 
ods m use a greedy heuristic algorithm. To test the 
stability of this method, we run this algorithm 100 times 
on the same force network and find that the fluctuation 
of maximal value of the modularity Q is within 1 %. In 
addition, we observe that the 15 largest communities 
consist of the same core group of particles: 70% of the 
same subset of particles are included in 90% of the iter¬ 
ations. We additionally And that fluctuations in Q are 
not accompanied by fluctuations in the morphology of 
the detected communities (to be quantified in ^4.2). 
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Fig. 2 Five example communities (7 = 3) and their convex 
hulls, for a packing at with /i = 0.3 and P = 10“^. Only 
communities containing more than 10 particles are shown. 

The choice of resolution parameter 7 controls the 
total number of communities identified, and also their 
morphology. For 7 < 1, optimizing Q favors large com¬ 
munities, while for 7 > 1 small communities dominate. 
To select the optimal value of 7 , we seek a figure of merit 
which quantifies the extent to which the detected com¬ 
munities take on a chain-like character: branched and 
sparse. We found that the technique used by Bassett 
et al. [T] for 2 D packings was ineffective in 3D systems. 
We therefore define a new figure of merit, the normal¬ 
ized convex hull ratio 


where Vp is the total volume of particles in the com¬ 
munity and Vhuii is the volume of the convex hull of 
the community. Branched/sparse communities will have 
lower values of He. To calculate Vhuiu we discretize each 
sphere as a 7 x 7 matrix of points and determine the 
convex hull using the Matlab boundary function. Fig.[^ 
shows example convex hulls. 

Fig.j^shows two examples of intermediate-size com¬ 
munities, shown in isolation to make them more visible. 
Note the chain-like structures dominating the commu¬ 
nities, providing a sparse structure with a low hull ratio. 
The interstices of such communities can be filled either 
by smaller communities, or by intercollated communi¬ 
ties which are also branched. 

For each packing, we calculate the mean hull ratio H 
by averaging the measured He weighted by the number 
of particles in each community, excluding communities 
which contain only one particle. To determine the opti¬ 
mal value of 7 to use in our analysis, we measure how H 
changes as a function of 7 across a range of pressures. 



Fig. 3 Examples of two medium size communities. Resolu¬ 
tion parameter 7 = 1.1. fi = 0.3. P = 10“^. 


As shown in Fig. [^, there is a clear minimum value 
of H which is approximately consistent across different 
values of P. In the analysis below, we utilize 7 = 1.1 in 
all cases. 

The effect of changing 7 can be understood by ex¬ 
amining the contribution each community makes to the 
modularity Q. The network force 

^C= 0) 

i,jec 

is the contribution to Q (Eq. from only the particles 
located in a particular community C. Its value increases 
due to both the normal forces in the community being 
large, and from the size (number of nodes) Se in the 
community. 

In Fig. [^, the individual communities are colored 
by their particular network force cTc- For small 7 , the 
Wij term dominates the sum in Eq. If 7 is small 
enough that Wij — jPij is always positive, then the 
largest value of Q is obtained by putting many particles 
in the same community {S{ci^ Cj) is nearly always 1). Eor 
larger 7 , the null model Pij will have more influence on 
the chosen communities, and the particular geometry 
and interparticle forces matter. If 7 is large enough that 
Wij —"fPij is always negative, then the optimal Q is zero 
by letting all (5 (q, Cj) be zero. In that case, the optimum 
value of Q is obtained when each community contains 
only a single particle. 

Eor the special case where 7 = 1, contact forces be¬ 
tween particles (Wij) are directly compared with the 
average contact force (Pij) in the system. The modu¬ 
larity Q is increased when more particles with multiple 
contact forces greater than / are included in commu¬ 
nities (force chains in our sense). The choice of 7 = 1 
is similar to finding force chains by thresholding at a 
minimum force, often set to be /. However, in contrast 
with thresholding methods, the modularity maximiza¬ 
tion method is flexible rather than binary. 
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Fig. 4 (a) Average hull ratio if as a function of resolution parameter 7 for various pressures, with /j,=0.3. (b) Community 
detection results from selected pressures and 7 (0.4, 0.8, 1, 1.2, 3). The color scale for each packing ranges from zero (deep 
blue) to the its maximum value of ac (red); many blue particles are of similar color but do not below to the same community. 
For clarity, communities with only 1 particle are hidden. 


4 Results 


Using these community detection methods, we describe 
how the force chain network changes as a result of 
both interparticle friction /i and confining pressure P. 
For each community, we consider three properties: the 
community size Sc, the network force ac (community 
strength), and the hull ratio He (community morphol¬ 
ogy). In all cases, community-detection is performed at 


fixed resolution parameter 7 = 1 . 1 , chosen as a com¬ 
promise value for the whole parameter regime. 

4.1 Community Size and Strength 

To illustrate the methods, we first examine the set of 20 
configurations with P = 10“^. In Fig. sample com¬ 
munity assignments are shown for all seven ji values. At 
low values of ji (top row), small communities dominate, 
while at large values of /i, there is typically a single 

















6 


Yuming Huang, Karen E. Daniels 



Fig. 5 Community assignments at all 7 values of interparticle 
friction /j, and P = 10“^. Color indicates the network force ac 
of each community. The color scale for each packing ranges 
from zero (deep blue) to the its maximum value of ac (red). 


large community near the top and many smaller and 
weaker communities at the bottom. (The multiple low- 
cTc communities all have similar values of ac and thereby 
appear to be the same community (by color) although 
they are not.) This observation is similar to prior work 
on the effect of friction coefficient fi on jamming proper¬ 
ties of packings [32] in which the bulk packing fraction 
and coordination number gradually decrease as /i in¬ 
creases from 0 and they saturate when /i is larger than 
1. This saturation is also reflected in the cumulative 
distribution figures we examine below. 

As shown in Fig. we observe that Sc and ac obey 
an approximately linear relationship. We examine their 
relationship under all p and P settings and find the 
same approximate relationship. (Since the mean pres¬ 
sure was already normalized in writing the weighted 
adjacency matrix W, we do not expect a trend in the 
magnitude of ac-) 

4-1.1 Friction-Dependence 



Fig. 6 Scatter plot of community size Sc and network force 
(7c, with each point representing a single community. Data is 
from a single simulation with ^ = 0.3 and P = 10“^. 



Fig. 7 Comparison of the cumulative distribution of commu¬ 
nity size Sc as a function of ii, where each plot represents the 
average over 20 simulations at each P. 
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Fig. 8 Comparison of the cumulative distribution of network 
force (7c as a function of /x, where each plot represents the 
average over 20 simulations at each P. 


To understand the friction-dependence, we consider the 
cumulative distribution function (CDF) of both Sc and 
network force ac as a function of p at fixed P. As shown 
in Figs. and both quantities show similar behav¬ 
ior, as expected given the strong correlation show in 
Fig-H For large p, we observe an approximately expo¬ 
nential distribution. Remarkably, the steepness of the 
distribution as a function of p has opposite trends at 


low and high pressure: For P < 10^, the CDF steepens 
as p increases (fewer large/strong communities), while 
for P > 10^, the CDF instead steepens as p decreases. 
Thus, P* ~ 10^ represents a crossover value between 
two distinct behaviors. Below, we will explore how the 
heterogeneity of forces (shown illustratively in Fig. 
causes this effect. 
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Fig. 9 Comparison of the cumulative distribution of commu¬ 
nity size Sc as a function of P, where each plot represents the 
average over 20 simulations at each fi. 


Fig. 10 Comparison of the cumulative distribution of net¬ 
work force ac as a function of P, where each plot represents 
the average over 20 simulations at each /j. 


4 . 1.2 Pressure-Dependence 

Fig.i and show the same CDF data, rearranged 
to highlight P-dependence at fixed p. This configura¬ 
tion highlights the existence of a low-friction regime 
distinct from the frictional regime, with a transition 
near fp ^ 0.1. For /i < /i*, the CDFs become much 
steeper as P is increased. This indicates that the sys¬ 
tem’s forces are becoming more homogeneous at high 
pressure, as expected [33[|34]. In contrast, simulations 
performed at /i > /i* (the frictional regime) show only 
weak pressure-dependence, with the large-dc tails fluc¬ 
tuating. This may be due either to insufficient statistics, 
or to changes in the heterogeneity of the system, to be 
discussed in the next section. 

These CDFs of Sc and ctc are similar to those ob¬ 
served in a previous study of force chains in 2D sys¬ 
tems [T] using a similar community-detection technique. 
There, the community size distribution was also ex¬ 
ponential, and here we found that the network force 
was exponential as well. In addition, both studies saw 
that communities are more compact at high pressure. 
This pressure-dependence is in contrast to the work of 
Navakas et al. m, in which it was observed that com¬ 
munity size increases as pressure increases. A key dis¬ 
tinction between the two studies is the choice of null 
model: they used the standard Newman-Girvan null 
model [35l [3T] rather than a geometric null model [T] , 
resulting in domain-like communities. 


4 . 1.3 Network Homogeneity 

We have observed that there is a crossover in commu¬ 
nity size and strength for both pressure (P* ~ 10^) 
and friction (/i* « 0.1), and that this effect appears 
to be connected to the homogeneity of the force net¬ 
work. To examine this in more detail, we consider the 
vertical gradient in the community size Sc and its re¬ 
lationship to the relative importance of horizontal vs. 
vertical forces. 

Fig. shows the spatial distribution of average 
community size {Sc) as a function of the vertical po¬ 
sition ^ within the sample, for each pair of (/i, P) pa¬ 
rameters. Averages are calculated on the particle-scale: 
within horizontal slice of thickness d, we average the 
Sc of all particles whose centers are within that bin. 
We observe that the plots fall into three distinct types: 
negative slope (colored red, largest communities at the 
bottom), an almost vertical distribution (colored yel¬ 
low, community size evenly distributed), and positive 
slope (colored blue, largest communities at the top). 
As expected from Fig. the corresponding plot for cjc 
is very similar (not shown). 

Note that the most homogeneous communities ap¬ 
proximately correspond to the P* = 10“^ crossover vis¬ 
ible in Fig. suggesting that spatial gradients are im¬ 
portant. For /i > /i*, the largest pressures used were 
able to reverse the gradient, moving the largest gradi¬ 
ents from the bottom to the top of the packing. Note 
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Fig. 11 Vertical distribution of average community size (Sc) 
at all values of fi and P settings, averaged over all 20 simula¬ 
tions. All axis scales are the same. Background colors indicate 
the slope of the plot: red for negative slope, blue for positive 
slope, and yellow for uniform community size. 
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Fig. 12 Vertical distribution of average (vertical, horizontal 
components of the) normal force f at all values of fi and P 
settings, averaged over all 20 simulations. All axis scales are 
the same. 


that these two kinds of gradients distinguish the similar- 
width distributions at P = 0 and P = 10“^ in Fig. 
when fji > fp. This non-monotonic dependence of het¬ 
erogeneity on pressure was unexpected. 

To understand how this effect arises, we consider the 
relative importance of horizontal and vertical normal 
interparticle forces as a function of 2 ;. As done for (Pc), 
we calculate the particle-scale average of the horizon¬ 
tal and vertical components of the vector normal force 
f. As shown in Fig. the interparticle forces mostly 
increase with depth, as would be expected for gravita¬ 
tional loading. (Due to the periodic boundary condi¬ 
tions in the lateral direction, the forces cannot saturate 
due to the Janssen effect [36].) This gravitation-loading 


10^ 


10^ 

10^ 


0 0.2 0.4 0.6 0.8 1 1.2 1.4 

Fig. 13 Scatter plot of hull ratio He as a function of commu¬ 
nity size Sc for a single simulation at P = 10“^ and /i = 0.3. 
Values of ife > 1 are possible because we approximate spheres 
as polyhedra in finding the convex hull. 
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regime approximately corresponds with the red-shaded 
plots in Fig. [U and is also visible in the middle column 
of Fig. Ef, where a big, strong community forms at the 
bottom part of the packing. In contrast, for P > P* the 
forces are more spatially homogeneous; similar effects 
have been seen by Makse et al. [aSHaU- For high P and 
/i (the blue-shaded plots), the vertical forces first be¬ 
come more uniform with depth, but eventually develop 
a force-excess at the top of the packing. This high-force 
region corresponds to the large communities shown in 
the bottom row of Fig. Between these two extreme 
cases, there is a regime in which the community-size 
distributions are quite homogeneous (the yellow-shaded 
plots in Fig. 11). This regime does not precisely corre¬ 
spond with the most homogeneous force distributions 
shown in Fig. 12, suggesting that community-detection 
is sensitive to small changes in the interparticle forces. 


4.2 Community Morphology 


While visual inspection of force chain morphology is 
possible in 2D systems, it is harder to observe such 
changes within 3D systems (see Fig. [^. Therefore, a 
key benefit of using the geographical null model (Eq. 
to detect the communities of particles which form the 
backbone of the system is to provide a way to quantify 
changes in the force chain network. The hull ratio He 
(Eq.§ measures the degree to which the communities 
are sparse/branched. In this section, we characterize 
how He changes as a function of p and P. As shown 
in Fig.[^ we observe that the largest communities are 
also the most branched (low He). An exception to this 
trend occurs when a large, strong community forms at 


the top of packing (the blue-shaded plots in Fig. 11). 

Fig. shows the cumulative distributions of hull 
ratio Pc: organized by pressure. Because > 90% of the 
communities contain only a single particle, these com- 
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Fig. 14 Comparison of the cumulative distribution of hull 
ratio He as a function of fi, where each plot represents the 
average over 20 simulations at different P. Low He values 
correspond to branched and sparse communities, while high 
He values correspond to dense communities. 


tion (/i* 0.1). In addition, this effect appears to be 

connected to the homogeneity of the force network. 

It is our hope that this technique will prove useful 
for investigating the statistical properties of force chain 
networks, by identifying the most important communi¬ 
ties of particles. While we have not included tangential 
forces in this study, including will likely be necessary 
for addressing questions of mechanical stability. 
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munities are not shown on the plots and have been ex¬ 
cluded from calculations of average hull ratio (including 
in Fig. [^. For P > P*, we observe that the cumula¬ 
tive distributions are sensitive to p; below P*, they are 
/i-independent. In the /i-dependent regime, we see that 
larger frictional forces contribute to finding more chain- 
like communities (low He). Conversely, it is also true 
that for high /i, the CDF of is more sensitive to P. 
This is consistent with studies in two dimensions [1], 
where the community shape for a frictionless packing 
was less sensitive to pressure than in frictional pack¬ 
ings. 

5 Conclusions 

In this paper, we have shown that community detection 
methods can be successfully applied to 3D granular ma¬ 
terials. We define a new quantity, the hull ratio, which 
characterizes the degree of branching within a commu¬ 
nity. This quantity allows us to optimize the community 
detection process by identifying a resolution (7 = 1 . 1 ) 
where we detect the most-branched features of the sys¬ 
tem. This resolution is in approximate agreement with 
observations in 2 D granular systems [T], and is sensi¬ 
ble given the normalization of the weighted adjacency 
matrix W. 

For packings generated over a range of interparticle 
friction fi and pressure P, we characterize the detected 
communities in terms their size, strength, and hull ra¬ 
tio. The first two are found to be largely redundant, 
and all three depend on /i and P. We find that, as in 
2D systems [1], the size and strength exhibit approxi¬ 
mately exponential distributions. Using these measures, 
we observe that there is a crossover in community size 
and strength for both the pressure (P* 10^) and fric- 
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